##### ADDITIONAL ANALYSES #########----------------------------------------------

# Script to generate all other additional analyses including all robustness tests.

# Session info
# R version 4.1.0 (2021-05-18)
# Platform: x86_64-apple-darwin17.0 (64-bit)
# Running under: macOS Big Sur 11.6.8

#-------------------------------Preparation-------------------------------------

# Load data
load("Data/vdem_panel_final.Rda")

# prepare data e.g. log and lags
vdem_panel_final <- vdem_panel_final %>% 
  mutate(oil_gas_log = log(oil_gas_value_nom + 0.01),
         v2x_polyarchy_lag = lag(v2x_polyarchy),
         e_polity2_lag = lag(e_polity2)) %>% 
  arrange(country_name, year) %>% 
  group_by(country_name) %>% 
  mutate(across(c(v2x_polyarchy, v2cademmob, v2caautmob,
                  gdppc_log_lag, growth_rate_3yrs_lag, oil_gas_log,
                  pop_log_lag, e_peaveduc, ethnic_frac,
                  intrastate_conflict, v2x_ex_party,
                  v2x_ex_military, v2x_ex_hereditary), ~ .x - lag(.x),
                .names = "{.col}_diff"),
         year_fac = factor(year)) %>% 
  ungroup()

#-------------Polity as DV (Table B2)-------------------------------------------

# pooled model
mass.pooled.pol <- fixest::feols(e_polity2 ~ v2cademmob_osp_lag + v2caautmob_osp_lag +
                                   gdppc_log_lag + growth_rate_3yrs_lag +
                                   oil_gas_log + pop_log_lag + 
                                   e_peaveduc_lag + 
                                   ethnic_frac + intrastate_conflict +
                                   v2x_ex_party_lag +  v2x_ex_military_lag + 
                                   v2x_ex_hereditary_lag,
                                 data = vdem_panel_final,
                                 panel.id = c("country_name", "year"),
                                 vcov = "DK")

#summary(mass.pooled.pol)

# TW FEs 
mass.fixed.pol <- fixest::feols(e_polity2 ~ v2cademmob_osp_lag + v2caautmob_osp_lag +
                                  gdppc_log_lag + growth_rate_3yrs_lag +
                                  oil_gas_log + pop_log_lag + 
                                  e_peaveduc_lag + 
                                  ethnic_frac + intrastate_conflict +
                                  v2x_ex_party_lag +  v2x_ex_military_lag + 
                                  v2x_ex_hereditary_lag |
                                  country_name + year_fac,
                                data = vdem_panel_final,
                                panel.id = c("country_name", "year"),
                                vcov = "DK")

#summary(mass.fixed.pol)

# LDV model 
mass.LDV.pol <- fixest::feols(e_polity2 ~ e_polity2_lag + v2cademmob_osp_lag + v2caautmob_osp_lag +
                                gdppc_log_lag + growth_rate_3yrs_lag +
                                oil_gas_log + pop_log_lag + 
                                e_peaveduc_lag + 
                                ethnic_frac + intrastate_conflict +
                                v2x_ex_party_lag +  v2x_ex_military_lag + 
                                v2x_ex_hereditary_lag |
                                country_name + year,
                              data = vdem_panel_final,
                              panel.id = c("country_name", "year"),
                              vcov = "DK")

#summary(mass.LDV.pol)

# FD model 
mass.FD.pol <- plm::plm(e_polity2 ~ + v2cademmob_osp_lag + v2caautmob_osp_lag +
                          gdppc_log_lag + growth_rate_3yrs_lag +
                          oil_gas_log + pop_log_lag + 
                          e_peaveduc_lag + 
                          ethnic_frac + intrastate_conflict +
                          v2x_ex_party_lag +  v2x_ex_military_lag + 
                          v2x_ex_hereditary_lag,
                        index = c("country_name", "year"),
                        model = "fd", data = vdem_panel_final)


# Export regression results to latex
texreg::htmlreg(file = "Output/table_B2.html", list(mass.pooled.pol, mass.fixed.pol,
                    mass.LDV.pol, mass.FD.pol),
               booktabs = T,
               omit.coef = "(Intercept)", 
               custom.coef.names = c("Mobilization for democracy",
                                     "Mobilization for Autocracy",
                                     "GDP per capita (log)",
                                     "GDP growth rate",
                                     "Gas and oil production (log)",
                                     "Population (log)",
                                     "Avg. years of education",
                                     "Ethnic fractionalization",
                                     "Intrastate conflict",
                                     "Ruling party dimension",
                                     "Military dimension",
                                     "Hereditary dimension",
                                     "Democracy (lag)"))

#-------------Extended controls (Table B3)--------------------------------------

# pooled model
mass.pooled_ext <- fixest::feols(v2x_polyarchy ~ v2cademmob_osp_lag + v2caautmob_osp_lag +
                                   gdppc_log_lag + growth_rate_3yrs_lag +
                                   oil_gas_log + pop_log_lag + 
                                   e_peaveduc_lag + 
                                   ethnic_frac + intrastate_conflict +
                                   v2x_ex_party_lag +  v2x_ex_military_lag + 
                                   v2x_ex_hereditary_lag +
                                   election_year +
                                   coup_attempt +
                                   theta_mean,
                                 data = vdem_panel_final,
                                 panel.id = c("country_name", "year"),
                                 vcov = "DK")

#summary(mass.pooled_ext)

# TW FEs 
mass.fixed_ext <- fixest::feols(v2x_polyarchy ~ v2cademmob_osp_lag + v2caautmob_osp_lag +
                                  gdppc_log_lag + growth_rate_3yrs_lag +
                                  oil_gas_log + pop_log_lag + 
                                  e_peaveduc_lag + 
                                  ethnic_frac + intrastate_conflict +
                                  v2x_ex_party_lag +  v2x_ex_military_lag + 
                                  v2x_ex_hereditary_lag +
                                  election_year +
                                  coup_attempt +
                                  theta_mean|
                                  country_name + year_fac,
                                data = vdem_panel_final,
                                panel.id = c("country_name", "year"),
                                vcov = "DK")

#summary(mass.fixed_ext)

# LDV model 
mass.LDV_ext <- fixest::feols(v2x_polyarchy ~ v2x_polyarchy_lag + v2cademmob_osp_lag + v2caautmob_osp_lag +
                                gdppc_log_lag + growth_rate_3yrs_lag +
                                oil_gas_log + pop_log_lag + 
                                e_peaveduc_lag + 
                                ethnic_frac + intrastate_conflict +
                                v2x_ex_party_lag +  v2x_ex_military_lag + 
                                v2x_ex_hereditary_lag+
                                election_year +
                                coup_attempt +
                                theta_mean |
                                country_name + year,
                              data = vdem_panel_final,
                              panel.id = c("country_name", "year"),
                              vcov = "DK")

#summary(mass.LDV_ext)

# FD model 
mass.FD_ext <- plm::plm(v2x_polyarchy ~ + v2cademmob_osp_lag + v2caautmob_osp_lag +
                          gdppc_log_lag + growth_rate_3yrs_lag +
                          oil_gas_log + pop_log_lag + 
                          e_peaveduc_lag + 
                          ethnic_frac + intrastate_conflict +
                          v2x_ex_party_lag +  v2x_ex_military_lag + 
                          v2x_ex_hereditary_lag +
                          election_year +
                          coup_attempt +
                          theta_mean, 
                        index = c("country_name", "year"),
                        model = "fd", data = vdem_panel_final)


# Export regression results to latex
texreg::htmlreg(file = "Output/table_B3.html", list(mass.pooled_ext, mass.fixed_ext, mass.LDV_ext, mass.FD_ext),
               booktabs = T,
               omit.coef = "(Intercept)", 
               custom.coef.names = c("Mobilization for democracy",
                                     "Mobilization for autocracy",
                                     "GDP per capita (log)",
                                     "GDP growth rate",
                                     "Gas and oil production (log)",
                                     "Population (log)",
                                     "Avg. years of education",
                                     "Ethnic fractionalization",
                                     "Intrastate conflict",
                                     "Ruling party dimension",
                                     "Military dimension",
                                     "Hereditary dimension",
                                     "Elections",
                                     "Coup attempt",
                                     "Human rights protection",
                                     "Democracy (lag)"))

#-------------EDI components (Table B4)--------------------------------------

# Freedom of Expression and Alternative Sources of Information 
fixed_freeexp <- fixest::feols(v2x_freexp_altinf ~ v2cademmob_osp_lag + v2caautmob_osp_lag +
                                 gdppc_log_lag + growth_rate_3yrs_lag +
                                 oil_gas_log + pop_log_lag + 
                                 e_peaveduc_lag + 
                                 ethnic_frac + intrastate_conflict +
                                 v2x_ex_party_lag +  v2x_ex_military_lag + 
                                 v2x_ex_hereditary_lag |
                                 country_name + year_fac,
                               data = vdem_panel_final,
                               panel.id = c("country_name", "year"),
                               vcov = "DK")

#summary(fixed_freeexp)

# Freedom of Association
fixed_freeassoc <- fixest::feols(v2x_frassoc_thick ~ v2cademmob_osp_lag + v2caautmob_osp_lag +
                                   gdppc_log_lag + growth_rate_3yrs_lag +
                                   oil_gas_log + pop_log_lag + 
                                   e_peaveduc_lag + 
                                   ethnic_frac + intrastate_conflict +
                                   v2x_ex_party_lag +  v2x_ex_military_lag + 
                                   v2x_ex_hereditary_lag |
                                   country_name + year_fac,
                                 data = vdem_panel_final,
                                 panel.id = c("country_name", "year"),
                                 vcov = "DK")

#summary(fixed_freeexp)

# Suffrage
fixed_suffr <- fixest::feols(v2x_suffr ~ v2cademmob_osp_lag + v2caautmob_osp_lag +
                               gdppc_log_lag + growth_rate_3yrs_lag +
                               oil_gas_log + pop_log_lag + 
                               e_peaveduc_lag + 
                               ethnic_frac + intrastate_conflict +
                               v2x_ex_party_lag +  v2x_ex_military_lag + 
                               v2x_ex_hereditary_lag |
                               country_name + year_fac,
                             data = vdem_panel_final,
                             panel.id = c("country_name", "year"),
                             vcov = "DK")

#summary(fixed_suffr)

# free and fair elections
fixed_freefair <- fixest::feols(v2xel_frefair ~ v2cademmob_osp_lag + v2caautmob_osp_lag +
                                  gdppc_log_lag + growth_rate_3yrs_lag +
                                  oil_gas_log + pop_log_lag + 
                                  e_peaveduc_lag + 
                                  ethnic_frac + intrastate_conflict +
                                  v2x_ex_party_lag +  v2x_ex_military_lag + 
                                  v2x_ex_hereditary_lag |
                                  country_name + year_fac,
                                data = vdem_panel_final,
                                panel.id = c("country_name", "year"),
                                vcov = "DK")

#summary(fixed_freefair)

# elected officials index
fixed_elecoff <- fixest::feols(v2x_elecoff ~ v2cademmob_osp_lag + v2caautmob_osp_lag +
                                 gdppc_log_lag + growth_rate_3yrs_lag +
                                 oil_gas_log + pop_log_lag + 
                                 e_peaveduc_lag + 
                                 ethnic_frac + intrastate_conflict +
                                 v2x_ex_party_lag +  v2x_ex_military_lag + 
                                 v2x_ex_hereditary_lag |
                                 country_name + year_fac,
                               data = vdem_panel_final,
                               panel.id = c("country_name", "year"),
                               vcov = "DK")

texreg::htmlreg(file = "Output/table_B4.html", list(fixed_elecoff, fixed_freeassoc, fixed_freeexp, fixed_freefair,
            fixed_suffr), booktabs = T,
       omit.coef = "(Intercept)", 
       custom.coef.names = c("Mobilization for democracy",
                             "Mobilization for Autocracy",
                             "GDP per capita (log)",
                             "GDP growth rate",
                             "Gas and oil production (log)",
                             "Population (log)",
                             "Avg. years of education",
                             "Ethnic fractionalization",
                             "Intrastate conflict",
                             "Ruling party dimension",
                             "Military dimension",
                             "Hereditary dimension"))


#-------------Additional lagged DVs (Table B5)----------------------------------
# FD model 
mass.FD_lags <- plm::plm(v2x_polyarchy ~  v2cademmob_osp_lag + v2caautmob_osp_lag + 
                           lag(v2cademmob_osp_lag,1) + lag(v2caautmob_osp_lag,1) +
                           lag(v2x_polyarchy, 1) +
                           lag(v2x_polyarchy, 2) +
                           lag(v2x_polyarchy, 3) +
                           gdppc_log_lag + growth_rate_3yrs_lag +
                           oil_gas_log + pop_log_lag + 
                           e_peaveduc_lag + 
                           ethnic_frac + intrastate_conflict +
                           v2x_ex_party_lag +  v2x_ex_military_lag + 
                           v2x_ex_hereditary_lag,
                         index = c("country_name", "year"),
                         model = "fd", data = vdem_panel_final)

#summary(mass.FD_lags)

texreg::htmlreg(file = "Output/table_B5.html", coeftest(mass.FD_lags, vcov=vcovHC(mass.FD_lags,type="HC0")))


#-------------FE counterfactual estimator (Figures B9 and B10)------------------
# devtools::install_github('xuyiqing/fastplm')
# devtools::install_github('xuyiqing/fect')

# dichotomizing treatment based on ordinal variable
vdem_fect <- vdem_panel_final %>% 
  mutate(demmob_bin_lag = ifelse(v2cademmob_ord_lag > 2, 1, 0),
         autmob_bin_lag = ifelse(v2caautmob_ord_lag > 2, 1, 0)) %>% 
  drop_na(v2cademmob_osp_lag)

# Counterfactual estimator for pro-democracy mobilization
out.fect_dem <- fect(v2x_polyarchy ~ demmob_bin_lag + autmob_bin_lag +
                       gdppc_log_lag + growth_rate_3yrs_lag +
                       oil_gas_log + pop_log_lag + 
                       e_peaveduc_lag + 
                       ethnic_frac + intrastate_conflict +
                       v2x_ex_party_lag +  v2x_ex_military_lag + 
                       v2x_ex_hereditary_lag, data = vdem_fect,
                     index = c("country_name","year"), na.rm = T,
                     cl = "country_name",
                     force = "two-way", method = "both", CV = T, 
                     vartype = "bootstrap", alpha = .05,
                     se = TRUE, nboots = 500, parallel = TRUE,
                     criterion = "gmspe", seed = 180722)

# Placebo test for pro-democracy mobilization
out.fect_dem_plac <- fect(v2x_polyarchy ~ demmob_bin_lag + autmob_bin_lag +
                            gdppc_log_lag + growth_rate_3yrs_lag +
                            oil_gas_log + pop_log_lag + 
                            e_peaveduc_lag + 
                            ethnic_frac + intrastate_conflict +
                            v2x_ex_party_lag +  v2x_ex_military_lag + 
                            v2x_ex_hereditary_lag, data = vdem_fect,
                          index = c("country_name","year"), na.rm = T,
                          force = "two-way", method = "both", CV = F, 
                          vartype = "bootstrap",alpha = .05,
                          cl = "country_name",
                          se = TRUE, nboots = 500, parallel = TRUE,
                          criterion = "gmspe",  placeboTest = T,
                          placebo.period = c(-5, 0), seed = 180722)

# Counterfactual estimator for pro-autocracy mobilization
out.fect_aut <- fect(v2x_polyarchy ~ autmob_bin_lag + demmob_bin_lag +
                       gdppc_log_lag + growth_rate_3yrs_lag +
                       oil_gas_log + pop_log_lag + 
                       e_peaveduc_lag + 
                       ethnic_frac + intrastate_conflict +
                       v2x_ex_party_lag +  v2x_ex_military_lag + 
                       v2x_ex_hereditary_lag, data = vdem_fect,
                     index = c("country_name","year"), na.rm = T,
                     force = "two-way", method = "both",
                     CV = TRUE, alpha = .05, 
                     cl = "country_name",
                     se = TRUE, nboots = 500, parallel = TRUE,
                     criterion = "gmspe", seed = 180722)

# Placebo test for pro-autocracy mobilization
out.fect_aut_plac <- fect(v2x_polyarchy ~ autmob_bin_lag + demmob_bin_lag +
                            gdppc_log_lag + growth_rate_3yrs_lag +
                            oil_gas_log + pop_log_lag + 
                            e_peaveduc_lag + 
                            ethnic_frac + intrastate_conflict +
                            v2x_ex_party_lag +  v2x_ex_military_lag + 
                            v2x_ex_hereditary_lag, data = vdem_fect,
                          index = c("country_name","year"), na.rm = T,
                          force = "two-way", method = "both", CV = F, 
                          alpha = .05,
                          cl = "country_name",
                          se = TRUE, nboots = 500, parallel = TRUE,
                          criterion = "gmspe",  placeboTest = T,
                          placebo.period = c(-5, 0), seed = 180722)

# Effects for demmob
plot(out.fect_dem, ylab = "Effect of dem. mobilization on democracy", 
     cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8)

ggsave(filename =  "output/Figure_B9a.pdf",
       plot = last_plot(), width = 5, height = 6, dpi = "retina")

plot(out.fect_dem_plac, ylab = "Effect of dem. mobilization on democracy", 
     cex.text = 0.8, stats = c("placebo.p","equiv.p"))

ggsave(filename =  "output/Figure_B9b.pdf",
       plot = last_plot(), width = 5, height = 6, dpi = "retina")

plot(out.fect_dem, "Average prediction error", 
     type = "equiv", 
     cex.legend = 0.8, main = "Testing Pre-Trend (IFEct)", cex.text = 0.8)

ggsave(filename =  "output/Fig_B9c.pdf",
       plot = last_plot(), width = 5, height = 6, dpi = "retina")

# Effects for autmob
plot(out.fect_aut, ylab = "Effect of aut. mobilization on democracy (IFEct)", 
     cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8)

ggsave(filename =  "output/Fig_B10a.pdf",
       plot = last_plot(), width = 5, height = 6, dpi = "retina")

plot(out.fect_aut_plac, cex.text = 0.8,
     ylab = "Effect of aut. mobilization on democracy",
     stats = c("placebo.p","equiv.p"))

ggsave(filename =  "output/Fig_B10b.pdf",
       plot = last_plot(), width = 5, height = 6, dpi = "retina")

plot(out.fect_aut, type = "equiv", ylim = c(-.1,.1),
     ylab = "Testing Pre-Trend (IFEct)",
     cex.legend = 0.6, cex.text = 0.8)

ggsave(filename =  "output/Fig_B10c.pdf",
       plot = last_plot(), width = 5, height = 6, dpi = "retina")


#-------------EDI does not predict mobilization---------------------------------

edi.FD.dem <- plm::plm(v2cademmob ~ v2x_polyarchy_lag  + v2caautmob_osp_lag +
                         gdppc_log_lag + growth_rate_3yrs_lag +
                         oil_gas_log + pop_log_lag + 
                         e_peaveduc_lag + 
                         ethnic_frac + intrastate_conflict +
                         v2x_ex_party_lag +  v2x_ex_military_lag + 
                         v2x_ex_hereditary_lag,
                       index = c("country_name", "year"),
                       model = "fd", data = vdem_panel_final)

summary(edi.FD.dem)

edi.FD.aut <- plm::plm(v2caautmob ~ v2x_polyarchy_lag  + v2caautmob_osp_lag +
                         gdppc_log_lag + growth_rate_3yrs_lag +
                         oil_gas_log + pop_log_lag + 
                         e_peaveduc_lag + 
                         ethnic_frac + intrastate_conflict +
                         v2x_ex_party_lag +  v2x_ex_military_lag + 
                         v2x_ex_hereditary_lag,
                       index = c("country_name", "year"),
                       model = "fd", data = vdem_panel_final)

summary(edi.FD.aut)

#--------------Method of composition (Figure B8)--------------------------------

# load posteriors for autmob
post_autmob <- readRDS("Data/v2caautmob_z_stretched.rds")

# Z-sample matrix:
post_autmob <- as.data.frame(post_autmob)

# prepare for analysis
post_autmob_new <- post_autmob %>%  rownames_to_column() %>% 
  mutate(country_text_id = str_sub(rowname, 1, 3),
         year = as.numeric(str_sub(rowname, 5, 8)) +1) %>% 
  select(country_text_id, year, everything()) %>% 
  distinct(rowname, .keep_all = T) %>% 
  select(-rowname) %>% 
  rename_with(.fn = ~ str_replace(.x, "V", "autmob_"),
              .cols = starts_with("V")) %>% 
  drop_na()

# load posteriors for demmob
post_demmob <- readRDS("Data/v2cademmob_z_stretched.rds")

# Z-sample matrix:
post_demmob <- as.data.frame(post_demmob)

# Prepare for analysis
post_demmob_new <- post_demmob %>%  rownames_to_column() %>% 
  mutate(country_text_id = str_sub(rowname, 1, 3),
         year = as.numeric(str_sub(rowname, 5, 8)) + 1) %>% 
  select(country_text_id, year, everything()) %>% 
  distinct(rowname, .keep_all = T) %>% 
  select(-rowname) %>% 
  rename_with(.fn = ~ str_replace(.x, "V", "demmob_"),
              .cols = starts_with("V")) %>% 
  drop_na()

# merge with full dataset
vdem_full_posteriors1 <- left_join(vdem_panel_final, post_demmob_new, by = c("country_text_id", "year")) 
vdem_full_posteriors <- left_join(vdem_full_posteriors1, post_autmob_new, by = c("country_text_id", "year")) %>% 
  distinct(country_name, year, .keep_all = T) 

# Run separate models for all posteriors (pro-autocracy mobilization)
fe_reg_aut <- vdem_full_posteriors %>% 
  select(starts_with("autmob_")) %>% 
  map( ~  feols(v2x_polyarchy ~ v2cademmob_osp_lag + .x +
                  gdppc_log_lag + growth_rate_3yrs_lag +
                  oil_gas_log + pop_log_lag + 
                  e_peaveduc_lag + 
                  ethnic_frac + intrastate_conflict +
                  v2x_ex_party_lag +  v2x_ex_military_lag + 
                  v2x_ex_hereditary_lag |
                  country_name + year_fac,
                data = vdem_panel_final),
       panel.id = c("country_name", "year"),
       vcov = "DK") %>% 
  map( ~ broom::tidy(.x, conf.int = T, conf.level = 0.95))

# Run separate models for all posteriors (pro-democracy mobilization)
fe_reg_dem <- vdem_full_posteriors %>% 
  select(starts_with("demmob_")) %>% 
  map( ~  feols(v2x_polyarchy ~ v2caautmob_osp_lag + .x +
                  gdppc_log_lag + growth_rate_3yrs_lag +
                  oil_gas_log + pop_log_lag + 
                  e_peaveduc_lag + 
                  ethnic_frac + intrastate_conflict +
                  v2x_ex_party_lag +  v2x_ex_military_lag + 
                  v2x_ex_hereditary_lag |
                  country_name + year_fac,
                data = vdem_panel_final),
       panel.id = c("country_name", "year"),
       vcov = "DK") %>% 
  map( ~ broom::tidy(.x, conf.int = T, conf.level = 0.95))

# extract estimates from model lists

# pro-autocracy mobilization
estimates_automob <- map_dfr(fe_reg_aut, 2)[2,] %>% 
  pivot_longer(cols = autmob_1:autmob_1800)

# pro-democracy mobilization
estimates_demmob <- map_dfr(fe_reg_dem, 2)[2,] %>% 
  pivot_longer(cols = demmob_1:demmob_1800)

# Plot estimates
ggplot2::ggplot(data = estimates_automob, aes(x = value)) +
  geom_density() +
  xlab("Monte Carlo Estimates of Mobilization for Autocracy Coefficient") +
  ggtitle("Mobilization for Autocracy Kernel Density")

ggsave("output/Figure_B8b.pdf", width = 7, height = 5)

# Plot estimates
ggplot2::ggplot(data = estimates_demmob, aes(x = value)) +
  geom_density()+
  xlab("Monte Carlo Estimates of Mobilization for Democracy Coefficient") +
  ggtitle("Mobilization for Democracy Kernel Density")

ggsave("output/Figure_B8a.pdf", width = 7, height = 5)



